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ABSTRACT 

In  this  paper  the  scalability  of  a  computationally 
intensive,  multi-scale  theory  exhibiting  progressive 
damage  in  a  3-D  woven  composite  is  investigated.  The 
theory  is  based  on  evolving  some  fundamental  damage 
modes  in  a  representative  volume  element  (RVE)  of  the 
micro- structure  of  a  3-D  woven  composite.  The  evolving 
damage  modes  affect  the  local  stresses  in  the  composite’s 
actual  microstructure  and  eventually  the  overall  stresses 
in  the  composite.  This  situation  is  considered  in  the  RVE 
via  a  transformation  field  analysis  (TFA).  Since  the  model 
is  computationally  intensive,  it  benefited  by  having  its 
computations  made  parallel.  This  was  accomplished  by 
coupling  the  model  to  LS-DYNA3D,  an  explicit 
Lagrangian  finite  element  code.  The  coupled  RVE-TFA- 
LSDYNA3D  code  enabled  the  study  of  weave-level 
damage  progression  in  the  3-D  woven  composites  under 
general  impact  conditions  and  could  be  a  useful  design 
tool.  In  the  present  work,  the  scalability  of  the  coupling 
between  the  RVE-TFA  and  LS-DYNA3D  codes  is 
examined  and  results  are  presented. 

1.  INTRODUCTION 

Light  weight  armors  used  in  soldier  and  vehicle 
protection  applications  owe  their  superior  performance  to 
the  armor  configurations  which  are  carefully  selected 
from  combinations  of  material  constituents  (e.g.,  fibers, 
resins,  ceramics,  composites  etc.)  and  their  architectures 
(e.g.  woven  or  stitched  composites,  etc).  The  design 
process  aimed  at  arriving  at  these  configurations  seeks  to 
disperse  and  attenuate  the  blast/projectile  impact 
generated  high  amplitude  shock  waves  to  mitigate  the 
energy  and  momentum  transfer  into  soldiers  and  their 
vehicles.  The  need  for  increasing  the  impact  resistance 
requires  high  strength  materials  such  as  ceramics  and  the 
need  for  dissipating  impact  energy  requires  resins  and 
elastomers.  When  two  such  dissimilar  material  systems 
are  required  to  function  together,  they  require  yet  a  third 
type,  in  the  form  of  woven  composites  or  encapsulations, 


to  hold  them  together  giving  shape  and  rigidity  to  the 
eventual  armor  structure.  Under  impacts,  the  three 
constituent  material  systems  interact,  deform  and  fracture 
in  a  controlled  manner  while  also  undergoing  myriad 
damage  modes  depending  on  how  each  material 
constituent  experiences  the  overall  impact. 

Other  than  vast  experimental  data  sets,  no 
procedures  are  currently  available  to  designers  or  material 
scientists  that  tell  them  why  certain  combinations  of 
materials  and  geometry  are  better  at  providing  optimal 
energy  and  momentum  transfer,  or  how  the  incipient  local 
damage  modes  are  absorbing  the  impact.  Experiments 
reveal  the  dominant  damage  modes  in  the  armor 
constituents.  Designers  explain  them  in  traditional 
damage  mechanics  perspective  of  failures  on  principal 
stress  planes,  but  have  not  transferred  such  knowledge 
into  the  computational  tools  beyond  linear  displacement 
based  stress  resolution. 


Porting  the  damage  mechanics  understanding  of 
experimental  databases  into  explicit  time  integration 
based  non-linear  finite  element  modeling  (FEM)  has  not 
progressed  recently  owing  to  the  difficulties  in  modeling 
the  involved  damage  modes  across  the  length  scales  of  the 
material  systems.  By  using  averaged  properties  and  thus 
ignoring  length  scales,  traditional  computational 
approaches  often  focus  only  on  how  the  overall  impact 
response  is  affected  by  the  built-up  nature  of  an  armor 
configuration  but  not  on  how  the  local  damage  mechanics 
is  enabling  such  improvement.  What  is  needed  then 
within  the  FEM  is  a  method  to  explicitly  model  the 
individual  damage  modes,  i.e.  damage  mechanics  handled 
at  least  on  two  length  scales:  one  the  familiar  global  scale 
of  the  armor  and  the  other  the  level  of  local  micros- 
structure  of  armor  constituents. 


One  such  method  reported  by  Bahei-El-Din, 
Y.A.,  Rajendran,  A.  M.  and  Zikry,  M.  A.,  (2003)  uses  the 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

DEC  2008 

2.  REPORT  TYPE 

N/A 

3.  DATES  COVERED 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Scalable  Computing  Of  The  Mesh  Size  Effect  On  Modeling  Damage 

5b.  GRANT  NUMBER 

iucuiaiiiES  in  yyuvcn  mmui 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Comput.  and  Inform.  Sciences  Directorate,  US  Army  Research 

Laboratory  Aberdeen  Proving  Grounds,  MD  21005 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release,  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

See  also  ADM002187.  Proceedings  of  the  Army  Science  Conference  (26th)  Held  in  Orlando,  Florida  on  1-4 
December  2008,  The  original  document  contains  color  images. 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

uu 

18.  NUMBER 
OF  PAGES 

8 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


transformation  field  analysis  (TFA)  of  a  representative 
volume  element  (RVE)  of  a  3-D  woven  composite  to 
model  up  to  nineteen  different  modes  of  damage;  each 
mode  consisting  of  either  sliding,  splitting,  bulk,  or 
tension  failures  of  fiber,  matrix  and  inter-phase  materials. 
While  multi-scale  theories  have  been  used  before,  the 
emphasis  has  been  on  the  global  response,  i.e.  very  fine 
meshes  in  the  global  domain  and  very  coarse  meshes  in 
the  local  micro  domain.  Even  when  large  local  meshes  are 
used,  allowing  the  elements  to  totally  fracture  and  form 
debris  clouds  doesn’t  yield  clues  as  to  the  underlying 
damage  modes.  By  ignoring  material  fragmentation,  the 
RVE-TFA  method  takes  a  mechanics  based  view  of  the 
developing  damage  modes,  i.e.  a  damage  mode  once 
initiated  will  grow  unhindered  but  controls  and  modifies 
the  surrounding  stress  fields.  The  evolving  stress  field 
may  further  initiate  damage  elsewhere  in  the 
microstructure. 


The  earlier  focus  in  the  RVE-TFA  work  was  on 
obtaining  the  local  microscopic  stresses  and  relating  them 
to  overall  stress  increments.  The  multi-scale  nature  of 
RVE-TFA  enabled  RVE  meshes  to  be  separate  from  the 
global  meshes,  and  thus  lead  to  the  consideration  of  larger 
local  RVE  meshes  for  a  more  realistic  modeling  of  the 
local  composite  weave  details,  and  the  associated  stress 
and  strain  fields,  (Valisetty  et  al.,  2005).  But  because  the 
local  RVE-TFA  analysis  are  conducted  at  each  and  every 
global  time  step,  the  attendant  RVE  computations  become 
time  consuming  placing  a  premium  on  the  RVE  mesh 
size.  Although  some  convergence  was  shown  for  the 
overall  stresses  (Valisetty  et  al.,  2007),  the  role  of  the 
RVE  mesh  size,  i.e.  how  much  it  helps  to  have  a  larger 
mesh  vs.  a  smaller  one,  needed  to  be  examined  from  the 
perspective  of  how  the  local  RVE  stress  distributions  are 
affected  by  the  locally  spreading  micro-damage.  This 
issue  was  studied  in  detail  and  the  progression  of  damage 
modes  under  the  applied  overall  strain  increments  was 
reported  in  2006  ASC.  By  detailing  the  effect  of  local 
micro  mesh  size  on  the  overall  response  as  well  as  on 
evolutions  of  the  local  damage  modes,  the  numerical 
underpinning  for  the  RVE-TFA  was  computationally 
demonstrated  (Valisetty,  2008). 

In  the  present  work,  the  RVE-TFA  model  was 
introduced  into  LS-DYNA3D,  an  explicit  non-linear 
dynamic  finite  element  code,  as  a  user  defined  material 
subroutine.  The  ability  of  this  subroutine  to  model  the 
effect  of  the  progressions  of  a  select  number  of  damage 
modes  on  the  overall  stress-strain  behaviors  was 
demonstrated  in  a  report  (2006  ASC).  The  results 
demonstrating  the  scalable  nature  of  the  insertion  of  the 
RVE-TFA  model  into  the  LS-DYNA3D  finite  element 
code  are  presented  in  this  work. 


2.  RVE-TFA  COUPLING  TO  LS-DYNA3D 

RVE-TFA  is  added  as  a  user  defined  material 
subroutine  to  parallel  LS-DYNA3D.  The  computations  of 
the  global  mesh  are  handled  by  LS-DYNA3D  and  are 
spread  among  processors  using  global  domain 
decomposition.  At  integration  points  where  RVE-TFA 
code  is  used,  LS-DYNA3D  supplies  the  current  global 
strain  increments  to  the  RVE-TFA  user  defined  material 
subroutines,  and  in  turn  the  subroutines  perform  the  TFA 
analyses  and  supply  back  the  global  stress  increments  to 
LS-DYNA3D. 


The  coupling  of  RVE-TFA  to  LS-DYNA3D  is 
accomplished  such  that  computations  pertaining  to  the 
global  mesh  are  performed  by  LS-DYNA3D  remain 
parallel,  while  the  computations  in  the  RVE-TFA  code 
are  done  in  a  serial  manner  by  the  processors  on  which 
the  RVE-TFA  user  defined  material  sub-routines  are 
called.  The  data  required  by  the  RVE-TFA  user  material 
subroutines  are  read  by  the  processors  from  within  these 
subroutines  and  are  held  in  local  arrays.  The  processors 
take  turns  reading  the  data,  and  by  keeping  the  data  in 
local  arrays,  the  global  data  arrays  that  LS-DYNA3D  uses 
for  its  part  of  parallel  global  mesh  computations  are  not 
disturbed,  and  the  scalability  of  the  eventually  coupled 
RVE-TFA-LS-DYNA3D  computations  are  held  at  the 
level  of  the  LS-DYNA3D  global  computations.  After 
describing  the  RVE-TFA  inputs  and  the  nature  of  the  data 
that  this  code  requires  from  one  time  step  to  the  next  one, 
results  are  presented  to  demonstrate  the  scalability  of  the 
coupled  code. 


The  coupling  of  RVE-TFA  to  LS-DYNA3D  rests  on 
the  fact  that  an  RVE  can  be  used  to  compute  the  global 
stresses  as  sums  of  elastic  stresses  computed  directly  from 
the  global  strain  increments  that  are  known  at  the 
integration  points  of  a  global  mesh  and  some  self- 
equilibrating  local  transformation  fields  which  are 
significant  within  the  RVE.  Each  of  the  transformation 
fields  are  associated  with  the  progression  of  a  damage 
mode  within  the  RVE;  the  magnitude  of  each  one  of  such 
fields  are  determined  to  keep  the  local  stresses  at  levels  to 
maintain  orderly  progressions  of  the  damage  modes  as  per 
some  pre-selected  damage  criteria;  and  the  task  of 
keeping  the  local  RVE  stresses  at  progressive  damage 
levels  is  achieved  by  adjusting  the  local  stresses  to  restore 
local  equilibrium  of  the  RVE.  The  TFA,  or  the 
transformation  field  analysis,  refers  to  the  last  step  of 
local  RVE  stress  adjustment.  The  multi-scale  nature  of 
RVE-TFA  approach  is  self-evident  from  the  facts  that 
different  RVE’s  can  be  considered  simultaneously  at 
different  integration  points  of  a  global  mesh  and 
calculations  performed  in  one  RVE  are  independent  of  the 
calculations  performed  in  other  RVE’s. 


While  the  computations  that  occur  within  the  LS- 
DYNA3D  follow  the  well  known  explicit  time  integration 
schemes,  contact  algorithms,  etc.,  the  key  to  the  success 
of  the  RVE-TFA  coupling  to  LSDYNA-3D  is  the 
recognition  of  the  nature  of  the  computations  that  occur  in 
the  RVE-TFA  code,  the  data  required  for  conducting  such 
computations,  and  how  that  data  is  managed  and  updated 
between  the  time  integrations  that  occur  in  LS-DYNA3D 
code.  While  the  end  results  of  TFA  are  the  global  stress 
increments,  these  do  come  about  after  posing  and  solving 
the  TFA  problem  with  the  task  of  keeping  the  locally 
developing  damage  modes  under  progressive  control  from 
one  global  time  step  to  the  next  one.  Thus  there  is  a 
presumption  of  progressive  damage  growth.  This  aspect  is 
also  described  next  together  with  the  information  on  how 
the  data  necessary  for  tracking  these  damage  modes  is 
handled;  because  this  aspect  is  the  key  to  achieving  the 
scalable,  coupling  of  the  RVE-TFA  to  LS-DYNA3-D. 


2.1.  LOCAL  STRESS  EQUATIONS 

The  theoretical  development  of  RVE-TFA  is 
described  in  literature.  From  among  the  equations  of  the 
theory,  only  those  equations  are  presented  below  that  help 
the  description  of  the  data  that  is  necessary  for  computing 
the  local  stresses,  the  transformation  fields  that  correct  the 
local  stresses,  and  the  data  to  keep  track  of  the  growth  of 
the  local  damage  modes. 

Fig.  (1)  shows  a  typical  RVE  used  for  a  3-D  woven 
composite.  Overall  dimensions  of  the  RVE  are  based  on 
the  dimensions  apparent  in  the  weave’s  micrograph.  The 
number  of  sub-volumes  in  an  RVE  can  be  arbitrary  but  it 
is  convenient  to  have  each  RVE  sub-volume  belonging  to 
a  distinct  materials:  fiber  bundle,  matrix,  or  interface. 


Fig.  1.  A  Schematic  of  a  3-D  woven  composite  and 
corresponding  RVE 


Original  formulation  of  the  TFA  is  due  to  Dvorak 
(1992)  who  later  used  the  method  for  inelastic  composite 
materials  in  Dvorak  et  al.  (1994).  Dvorak  and  Zhang 
(2001)  used  the  TFA  to  analyze  damage  evolution  in  two- 
phase  composites.  Bahei-El-Din  et  al.  (2004)  extended 
the  TFA  method  to  the  study  of  plate  impact  problems  of 
3-D  woven  composites  by  considering  a  full  complement 


of  damage  modes  including  matrix  and  fiber  cracking, 
fiber  sliding,  interface  failure  and  debonding. 

A  brief  summary  of  the  relevant  equations 
(Bahei-El-Din  and  Boutros,  2003)  of  TFA  is  presented 
below  using  a  notation  of  boldface  lower  case  letters  for 
representing  the  (6x1)  stress/strain  vectors,  and  boldface 
upper  case  letters  for  the  corresponding 
stiffness/compliance  matrices.  The  global 
stresses/ strains,  a  and  e  ,  at  an  (integration)  point  where 
an  RVE  is  used,  are  weighted  volume  sums  of  local  sub¬ 
volume  stresses  and  strains.  The  local  sub-volumes  are 
assumed  to  undergo  elastic  deformation.  Deviations  from 
this  mode  of  deformation  are  treated  by  adding 
transformation  fields  to  account  for  the  progressing 
damage  modes.  A  part  of  the  local  sub-volume 
stresses/strains  can  be  related  to  the  global  applied  loads 
following  the  treatment  of  elastic  composite  aggregates 
by  Hill  (1963,  1965),  and  the  rest  to  a  superposition  of  all 
the  transformation  fields  originating  in  all  local  sub¬ 
volumes  within  the  RVE,  Dvorak  (1992)  as 

sr  =Ar  e+  IDA,ar  =Br  ct+  IFrsA,s  r=  1,...,Q  (1) 

s=l  8=1 

£r=MrCT,+  Pr  ,  <V=LrEr+Xr  (2) 

CT=Zcrar,  (3) 

r-1 

where  Q  is  the  number  of  sub-volumes  in  the  RVE,  Q  is 
the  number  of  ongoing  damage  modes  in  the  RVE,  cr  the 
sub-volume  fractions,  Ar  and  Br  the  stress/strain 
concentration  factors,  Drs  and  and  Frs  the  transformation 
influence  factors  for  the  rth  local  sub-volume,  Lrand  Mr 
elastic  stiffness  and  flexibility  matrices  and  pr  and 
X  r  =  -  L  r  p  r  are  transformation  strains  and  stresses. 

All  the  factors,  Ar ,  Br ,  Drs ,  Frs ,  Lr  and 
M  r  depend  upon  the  RVE  local  geometry  and  properties 
of  the  local  materials  and  are  determined  from  an  elastic 
analysis  of  the  RVE  mesh  and  are  supplied  as  input  to  the 
RVE-TFA-LS-DYNA3D.  The  stress  concentration  factors 
are  computed  as  statically  equivalent  to  the  global  stresses 
available  on  the  RVE.  For  example,  the  k'h  column, 
k=l,...,6,  of  the  stress  concentration  factor,  Br,  r  =  1, 
2,. .  ,,Q,  is  computed  by  giving  a  unit  value  to  the  kth  stress 
component  and  0  to  the  rest  of  the  6x1  overall  stress 

vector,  a .  This  is  done  in  the  usual  finite  element  sense 
by  assuming  a  linearly  varying  local  displacement  field  in 
the  RVE,  and  by  securing  the  RVE  against  rigid  body 
deformation,  (Dvorak  and  Teply,  1985). 


A  similar  procedure  is  also  used  for  computing  the 
transformation  stress  influence  factors,  Frs ,  by  realizing 
that  a  kth  column,  k=l,2,...,6,  of  Frs ,  r,s=l,2,...,Q,  is  a 
response  in  the  sub-volume  Vr  after  a  unit  stress,  A,k  =1, 
is  applied  in  the  sub-volume  Vs .  A  total  of  6Q  RVE  finite 
element  solutions  are  required  to  completely  evaluate  the 
transformation  factors.  Since  the  RVE  stiffness  matrix 
remains  the  same,  the  only  variants  are  the  effective  loads 
on  the  right  hand  side. 


Upon  entering  the  RVE-TFA  user  defined  material 
subroutines  for  the  first  time,  each  of  the  LS-DYNA3D 
computing  processors  read  the  factors,  Br  and  Frs ,  as 
input.  After  the  input  is  read,  the  only  data  that  the 
subroutines  keep  track  of,  from  one  global  time  step  to  the 
next  one,  is  the  data  that  help  them  to  compute  the 
transformation  fields.  Naturally  such  tracking  is  necessary 
to  bring  the  effect  of  local  material  damage  into  the  local 
stress  computations.  This  is  described  next. 


2.2.  COMPUTATION  OF  TRANSFORMATION 
STRESSES 


In  the  RVE-TFA  equations  (1-3)  presented  above,  the 
local  transformation  fields  are  used  as  corrections  to  the 
local  stresses  computed  from  the  global  strain  increments 
under  an  elastic  material  assumption.  When  the  RVE 
materials’  are  elastic,  there  is  no  need  for  these 
transformation  fields.  When  some  damage  develops  in 
one  or  more  of  RVE  sub-volumes,  the  transformation 
fields  are  used  to  offset  the  local  RVE  stresses  from  their 
elastic  values  to  restore  the  specific  damage  criteria,  in  the 
sense  of  radial  return  plasticity.  The  fact  that  stress  loss 
due  to  spreading  damage  can  be  handled  in  this  manner 
was  demonstrated  by  Bahei-El-Din  (1996),  Dvorak  and 
Zhang  (2001)  and  Bahei-El-Din  and  Botrous  (2003). 


The  transformation  stress  fields  are  found  by 
recasting  Eqs.  (1).  For  an  RVE  of  a  composite  which  is 
descritized  into  Q  sub-volumes,  Vr, r  =  1,2,..., Q,  the 
undamaged  local  elastic  stresses  are  found  from  the  stress 
concentration  factors,  Br ,  and  the  applied  overall  stress 
as  Br  ct  ,  the  first  parts  on  the  right  hand  side  of  Eq.  (1)2. 
Suppose  now  there  is  some  progressing  damage  in  one  or 
many  sub-volumes,  Vr,r  =  1,2,...,Q*,Q*  <  Q  of  the  RVE. 
The  computed  elastic  stresses  then  will  exceed  the 
underlying  strength  values  in  those  sub-volumes.  One  or 
several  damage  modes  may  be  progressing  in  those  sub¬ 
volumes.  Then  to  bring  the  computed  elastic  stresses  to 
conform  with  the  requirements  for  an  orderly  progression 
of  each  and  every  one  of  the  progressing  damage  modes, 


transformation  stress  fields,  equal  in  number  to  those 
causing  the  progressing  damage  modes,  are  introduced; 
these  are  the  second  parts  in  the  right  hand  side  of  Eq. 
(1)2, 


Essentially  Eqs.  (1)  are  statements  that  the  local 
stresses,  i.e.  the  sums  of  the  elastic  stresses  and  the 
transformation  stresses,  should  be  set  at  values  as  required 
by  the  damage  criteria  for  the  involved  damage  modes. 
This  may  result  in  the  development  of  damage  in  new 
sub-volumes  and  as  such  the  process  is  repeated  to 
compute  the  local  stresses  in  the  presence  of  all  possible 
damage  modes  caused  under  the  applied  overall  strain 
increment.  The  transformation  fields  are  evaluated  setting 
the  current  stresses  on  the  left  hand  side  of  Eqs.  (1)  at 
values  that  can  be  sustained  in  the  presence  of  progressing 
damage  modes. 


To  present  the  equations  and  later  the  damage  modes 
and  their  growth  and  criteria,  local  stresses  are  considered 
in  local  coordinate  systems  with  appropriate  co-ordinate 
transformation  matrices,  Rp .  The  local  coordinate 

systems  are  not  unique;  they  change  depending  on  sub¬ 
volume  material  orientations  and  the  specific  damage 
criteria.  Referring  to  such  local  coordinate  systems  and 

using  stress  ratios,  (pkp) ,  which  are  current  stresses 
divided  by  their  elastic  undamaged  values,  the  equations 
for  computing  transformation  stress  are  presented  as 
(Bahei-El-Din  et  ah,  2004): 


X  FMdiag(ak)Xn=-I-R-1diag«))Rp)BpCT,  p=l,2,„Q*, 
>1=1 


/ 1  if  0  <<()<1 
\  0  if  c|)  =  l  ’ 


(4) 

(5) 


where  diag(x)  is  a  (6  Q*  x6  Q* )  matrix.  The  stress  ratios, 
cpkp)  ,  k=l,2,...,6,  p  =  l,2,...,Q*,  Q* <Q,  imply  pre-  and 

post-damage  loading  as  follows:  cpkp)  =  1  for  a  stress 
component  that  has  yet  not  violated  the  onset  of  damage, 
cpkp)  =  0  indicates  complete  unloading,  and  0<(p<p'<l 
indicates  progressing  damage  with  partial  property  decay. 
If  tpkp)<l  for  all  k=l,2,...,6  and  for  all  Q*  sub-volumes, 
then  there  are  6  Q  equations  in  Eqs.  (4)  for  evaluating  the 
transformation  fields  Xr  which  number  to  6  Q  .  If  some  of 

the  local  stress  components  are  unaffected  by  the 
underlying  damage  criteria,  the  corresponding  stress 
ratios  would  equal  unity.  In  such  cases ak=0,  and  the 
corresponding  transformation  stresses  will  not  be  needed 
and  enough  equations  will  drop  out  from  Eqs.  (4).  The 


equations  available  after  this  consideration  are  solved  for 
the  transformation  stresses  that  are  remaining. 


While  the  6  Q*  x6  Q*  matrix  in  the  right  hand  side  of 
Eqs.  (5)  does  not  often  change,  the  left  hand  side  is 
updated  at  every  time  step  in  response  to  the  progression 
of  the  damage  modes.  As  described  next,  this  damage 
modes’  progression  is  iteratively  computed  with  the 
knowledge  of  local  stresses  at  the  beginning  of  the  current 
time  step  as  well  as  the  knowledge  of  how  these  local 
stresses  are  causing  materials  in  the  RVE  sub-volumes  to 
be  either  on  loading,  unloading  and  reloading  paths. 
Since,  in  coupling  the  RVE-TFA  code  to  LS-DYNA3D  as 
user  defined  material  subroutine,  only  the  overall  strain 
increments  and  stress  increments  are  passed  between  the 
codes,  the  data  necessary  in  the  form  of  stresses  from  the 
previous  time  step  and  material  loading/unloading  state 
are  written  to,  and  read  from,  the  external  files. 


2.3.  DAMAGE  MODES  AND  PROGRESSION 

Although  any  microstructure  can  be  considered 
with  the  aid  of  a  separate  local  RVE  mesh,  the  one 
considered  in  this  study  is  suitable  for  modeling  3-D 
woven  composites.  In  this  section,  the  damage  modes  that 
are  generic  to  this  type  of  composite  are  described  and 
particulars  are  presented  with  regards  to  their  numerical 
implementation  in  coupling  the  two  RVE-TFA  LS- 
DYNA3D  codes.  Described  next  are  the  numerical  data 
required  for  computing  the  damage  progression  and  its 
effect  on  local  stresses. 


2.3.1  Damage  modes 

The  RVE-TFA  code  (Bahei-El-Din  et  al,  2004) 
considered  damages  of  fiber  bundles,  interface  layers,  and 
matrix  but  specializes  the  progression  of  these  damages 
on  planes  specific  to  3-D  woven  fabrics.  Specific  damage 
modes  are  considered  in  each  material  system.  Four 
damage  modes  are  considered  in  the  fiber  bundles  as 
follows:  longitudinal  fiber  rupture,  transverse  inter-bundle 
fiber  sliding,  longitudinal  inter-bundle  fiber  sliding,  and 
transverse  inter-bundle  fiber  splitting.  Two  damage 
modes  are  considered  in  the  interface  layers  as  follows: 
shear  sliding  and  peeling.  The  bulk  matrix  damage  is 
assumed  to  be  terminal,  i.e.  without  any  recovery. 

Some  of  these  damage  modes  interact  with  each  other 
and  while  others  do  not.  For  example,  <j22,<j33,<j23  ,  the 
stresses  resolved  in  a  co-ordinate  frame  aligned  with  the 
axis  of  a  fiber  bundle,  affect  fiber  bundle  sliding  and 
splitting.  Therefore,  these  modes  are  deemed  interactive. 
On  the  other  hand,  the  two  interface  layer  damage  modes 
do  not  interact,  and  these  are  deemed  independent.  The 


stress  ratios,  <pk  ,  which  are  introduced  to  keep  track  of 
the  damage  progression,  become  unique  and  independent 
when  the  corresponding  damage  modes  do  not  interact; 
while  the  stress  ratios  involved  in  the  interacting  damage 
modes  are  not  unique.  For  such  interacting  modes,  the 
stress  ratio  which  controls  the  weakest  of  the  interacting 
damage  modes  is  selected. 


2.3.2  Material  strength  loss  modeling 

For  all  the  damage  modes  considered,  the  stress  and 
strain  allowable  curves  are  assumed  using  scalars 
quantities,  (,?,e),  that  can  be  computed  with  the  knowledge 
of  the  local  stresses  and  the  damage  modes.  These  curves 
are  initially  assumed  to  be  linear,  until  the  ultimate 
stresses  and  the  corresponding  elastic  strain  limits  are 
reached  (.S'u|t,  eel).  They  are  later  assumed  to  show 
softening.  The  softening  branches  are  assumed  to  take  a 
linear  or  a  nonlinear  form,  defined  symbolically  as  s=g(e- 
eei),  where  eei  <  e  <  euit,  and  eu\t  are  the  ultimate  strains. 
Unloading/reloading  follows  a  linear  path  between  the 
origin  and  the  (,?,e)  curves. 


2.3.3  Damage  Progression 

The  afore-mentioned  post- elastic  material  behavior  of 
the  RVE  sub-volumes  is  implemented  with  the  aid  of 
numerical  loading,  i.e  by  keeping  track  of  stresses  at  the 
beginning  of  the  current  time  steps,  and  how  the  current 
stresses  are  causing  loading,  unloading  or  re-loading  w.r.t 
all  the  active  damage  modes. 

For  a  given  strain  e(1),  at  the  ith  global  time  step,  the 
elastic  stress  is  computed  as  s,<l)=(s'uit/eei)e<1),  which  may  or 
may  not  fall  on  the  allowable  {s,e)  curve.  The  correct 

stress  magnitude,  and  hence  the  current  stress  ratio,  <p " 1 , 

is  determined  to  match  the  allowable  from  the  (y,e)  curve 
and  the  values  of  ,v(1),  e(1),  and(pll_1) ,  the  stress  ratio  at  the 
previous  step,  as  follows: 


undamaged  state  :  cp(I>  =1,  if  cp"  1)=1,  0<e(1)<eel 
damagereversal :  (pll,=  (p(1  ”,  if  cp(1_1)<l,  0  <e(1)<ed 
continued  damage : 


(p(l)  = 


g(e(l)  -eel)  ..  g(e(l)  -eel) 


JO 


if 


JO 


—  «P(1  Vej  -e<l)  <eult 


damagereversal : 


,(i)_m(i-i)  g(g(l>-gei)  "  m(i-D 


(pw=cp,w\if 


>(pvl_1\eei<e(i)<e, 


JO 


ult 


(6) 


(7) 


The  1st  condition  in  Eq.  (6)  represents  the  response  of  an 
undamaged  material,  while  the  1st  condition  in  Eq.  (7) 


describes  continued  loading  of  a  previously  damaged 
state.  The  remaining  conditions  above  represent  linear 
unloading/reloading  from  a  damaged  state.  Hence, 
damage  is  recorded  during  the  loading  history  through  the 
incremental  values  of  the  stress  ratio  cp£ 1 . 


The  stress  and  strain  pair  at  the  current  loading  step 
(s"}  has  different  interpretations  depending  on  the 

type  of  the  RVE  material  and  the  associated  damage 
mode.  In  what  follows,  these  quantities  are  related  to  the 
local  stress  and  strain  averages  for  a  fiber  bundle,  the  bulk 
matrix,  and  the  matrix  interface  layers.  It  is  assumed  that 

the  total  overall  stress  <j<0  at  the  current  loading  step  is 
known,  either  directly  by  following  a  defined,  stress- 
controlled  loading  path,  or  as  the  stress  Br  L  e"’ caused 
in  an  undamaged  RVE  by  a  strain-controlled  loading,. 
The  corresponding  local  stress  in  a  sub¬ 
volume  V  ,p=l, 2,..., Q,  of  the  RVE,  can  then  be  found 

from  Eq.  (1)2.  When  this  is  transformed  to  the  local 
coordinates  of  the  sub-volume,  and  the  result  is 
substituted  into  the  first  term  of  Eq.  (2)2,  the  local  stress, 

orr ,  and  strain,  T:n  ,  described  in  the  local  coordinate 
system  of  the  sub-volume,  are  found  at  loading  step  (i)  as 


The  composite  panel  was  made  of  a  3-D  woven  made 
of  a  S2  glass/ epoxy  system.  The  composite  was  evaluated 
to  have  the  micro-geometry  shown  in  Fig.  2. 


Fig.  3.  Cross-sections  through  (a)  the  warp  fiber  bundle 
and  (b)  through  the  weft  fiber  bundle. 


'  0) 


R  B  ct  V  s®  =Ma 


(i) 


r  =  l,2,...,Q  . 


(8) 


Matrix  Rr  relates  the  stresses  described  in  the  local  and 
the  overall  coordinate  systems,  and  Mr  is  the  elastic 
compliance  described  in  the  local  axes.  Expressions  for 
evaluating  s'"  in  terms  of  the  resolved  stress  d  and 
material  properties  eel  _  eult  _  slllt  for  the  above  listed 
damage  conditions  follow  composite  material  mechanics 
damage  concepts.  The  specific  ones  used  in  this  study  are 
from  Bahei-El-Din  et  al.,  2004. 


4.  RESULTS 

In  this  section,  numerical  results  are  presented  to 
show  the  scalability  of  the  coupling  of  the  RVE-TFA 
code  to  the  parallel  LS-DYNA3D.  For  this  purpose  a 
global  impact  problem  shown  in  Fig. 2  was  considered  in 
which  a  circular  glass  composite  panel  was  under  impact 
by  several  impactors,  all  traveling  at  154  m/s.  The  global 
impact  response  was  handled  by  the  parallel  LS- 
DYNA3D,  the  local  composite  response  was  computed  by 
the  RVE-TFA  code. 


As  seen  in  the  micrographs  of  Fig.  (3a),  the  fiber  bundles 
in  the  warp  and  weft  directions  have  a  rectangular  cross- 
section  with  aspect  ratio  of  about  4.0,  while  the  z-fiber 
bundle  has  a  square  cross-section.  For  the  composite,  a 
simple  RVE  idealization  with  an  approximate  material 
placement  of  the  fiber  bundles  in  the  warp,  weft  and  z 
directions,  and  the  resin  and  interfacial  layers  is 
considered  shown  in  Fig.  (1)2.  The  dimensions  of  the 
RVE  parallel  to  the  overall  in-plane  axes  were  taken  as 
12mm  and  12mm,  and  in  the  z-direction  as  2.2mm. 


For  fiber  bundles,  the  following  properties  are  used: 
El=59.7  GPa,  Ex=11.7  GPa,  vLT=0.248,  vTT=0.371, 
Glt=5  GPa,  Gtt=4.68  GPa,  cmax-L=3.36  GPa,  yT=0.0112, 
omax-T=0.08  GPa,  et=0.0068,  £max_T=0.0068,  xmax.L=0.057 
GPa,  xmax_T  =  0.048  GPa,  yL=0.01 14,  ymax.L=0.01 14,  xmax_ 
T=0.048GPa,  ymax_T=0.0112,  where  E,  G,  v,  a,  x, 
e  and  y  are  Young’s  modulus,  shear  modulus,  Poisson 
ratio,  normal  stress,  shear  stress,  normal  strain  and  shear 
strain,  respectively,  and  subscripts  L  and  T  indicate 
longitudinal  and  transverse  directions.  For  the  matrix,  the 
following  properties  are  used:  E=2.9  GPa,  =0.3, 
Omax=0.06  GPa,  Emax=0.021,  xmax=0.035  and  ymax  =0.0313. 


Having  thus  developed  the  RVE  of  Fig.  (1)2,  they  are 
used  for  the  present  study.  To  begin  the  study,  several 
RVE  meshes  were  considered  as  shown  in  Fig.  (4).  Once 


again  several  possibilities  existed  for  increasing  the  sub¬ 
volumes,  but  the  indicated  pre-ponderance  of  sub¬ 
volumes  in  the  z-direction  was  in  view  of  the  intended 
application  of  RVE-TFA  in  composite  armor  impact 
design  studies. 


150 

490 

294 

882 

Fig.  4.  Different  finite  element  meshes  for  the  RYE 


For  computing  the  composite’s  global  response,  the 
problem  was  modeled  with  a  global  mesh  of  5360  nodes 
and  4320  elements.  At  integration  points  of  this  global 
mesh,  LS-DYNA3D  called  the  coupled  RVE-TFA  code  to 
compute  the  local  stresses,  any  local  damage,  and  finally 
the  global  stress  increments.  Since  the  size  of  the  local 
RVE  mesh  can  influence  the  RVE-TFA  predictions,  and  a 
larger  mesh  can  improve  such  predictions,  the  tendency  is 
to  use  a  larger  RVE  mesh  for  computing  the  global  stress 
increments  at  integration  points.  But  increasing  the  local 
RVE  mesh  sizes  increases  not  only  the  numerical  cost  of 
the  local  computations  and  but  can  also  affect  the 
scalability  of  the  overall  solution.  A  recent  study 
(Valisetty  2007)  showed  that  a  RVE  mesh  size  of  490  is 
sufficient  to  capture  all  the  damage  modes  that  are 
important  in  the  3-D  woven  composites.  The  purpose  of 
the  present  study  was  then  to  show  that  RVE  mesh  sizes 
up  to  490  elements  can  be  used  for  TFA’s  in  a  global 
impact  problem  without  affecting  the  scalability  of  global 
computations. 

Table  1:  Mesh  size  vs.  size  of  the  set  of 

the  Transformation  Influence  Factors 


RVE  mesh 

No.  of 

Size  of  all  Frs 

Elements 

Nodes 

150 

252 

809712 

294 

448 

3111696 

490 

704 

8643024 

882 

1216 

28000656 

Among  the  input  that  the  TFA  analysis  requires,  i.e., 
the  Br  and  Frs  factors  in  Eq.  1,  the  input  of  Frs  is 
voluminous  as  shown  in  Table  1.  The  very  first  time  after 
entering  the  user  defined  RVE-TFA  material  subroutines, 
each  of  the  parallel  compute  processors  were  made  to 
read  this  data  one  after  other.  After  subtracting  the  wait 
times  for  reading  this  input  data,  the  global  solution  times 
for  rumiing  70  global  time  steps  were  noted  and  were 
plotted  in  Fig.  5.  The  multi-processor  run  times  for 
solutions  with  each  of  the  four  RVE  mesh  sizes  were  non- 
dimensionalzed  with  respect  to  the  corresponding  serial, 
i.e.  one  processor,  mn  times. 


Fig.  5  Scalability  results 

For  all  the  four  RVE  mesh  size,  there  is  good 
scalability  for  up  to  20  processors.  The  scalability  reduces 
for  nins  with  higher  number  of  processors,  but  the  drop 
from  ideal  scalability  is  not  much.  Even  this  drop  can  be 
attributed  to  the  fact  that  on  some  processors  local 
damage  was  not  an  issue  and  the  respective  processors 
were  finishing  their  local  computations  earlier. 
Considering  the  fact  that  before  this  coupling  to  the 
parallel  LS-DYNA3D,  the  effort  of  computing  the 
coupled  global-local  solutions  with  one  processor  runs 
were  taking  days  for  problems  with  larger  global  meshes 
and  local  RVE  meshes  of  size  490-882,  the  drop  in 
scalability  for  larger  number  of  processor  runs  is  not  a  big 
numerical  cost;  the  present  study  should  enable 
undertaking  of  investigations  into  salient  damage  modes 
and  their  growth  in  real  impact  situations  for  armor 
composites. 


5.  CONCLUSIONS  AND  FUTURE  WORK 

In  this  paper,  a  computationally  intensive  RVE-TFA 
code  used  to  study  the  progression  of  damage  modes  in  3- 
D  woven  composite  microstructures  was  coupled,  in  the 
manner  of  a  user  defined  material  subroutine  to  provide 
overall  stress  increments  to,  the  parallel  LS-DYNA3D  a 
Lagrangian  explicit  code  used  in  many  armor  composite 
impact  studies. 

The  coupling  of  RVE-TFA  to  LS-DYNA3D  was 
accomplished  such  that  the  computations  pertaining  to  the 
global  mesh  are  performed  by  LS-DYNA3D  and  remain 
parallel,  while  the  computations  of  the  RVE-TFA  code 
were  done  in  a  serial  manner  by  the  processors  on  which 
the  RVE-TFA  user  defined  material  sub-routines  are 
called.  With  the  help  of  a  brief  summary  of  equations,  the 
various  inputs  required  for  conducting  the  TFA  were 
identified  and  a  description  was  provided  as  to  how  these 
inputs  were  read  and  maintained  by  the  individual 
processors.  The  numerical  implementation  of  the  damage 
modes’  progression  was  described  together  with  the 
description  of  the  data  required  for  this  purpose  and  how 
such  data  was  handled  by  the  processors  from  one  global 
time  step  to  the  next. 

With  the  aid  of  a  global  impact  on  a  3-D  woven 
composite,  the  scalability  of  the  coupling  of  the  two  codes 
was  demonstrated.  The  results  of  this  study  should  enable 
the  undertaking  of  detailed  studies  of  damage  modes’ 
propagation  in  3-D  woven  composites  in  future. 
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